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ABSTRACT 


The thesis describes an analysis of a Magnetic Anomaly Detection (MAD) buoy 
field. The analysis was based on estimates of the probability a submarine would be 
detected during an encounter with a MAD buoy. The estimates were determined by 
using a simulation of both a crosscorrelation detection system and a square law 
detection system. Using these estimates, a random search model was used to 
determine a lower bound on the search effectiveness of a MAD buoy field. In doing 
this, the effects of false alarm rate, sweep width, noise level, buoy depth, and target 
displacement were analyzed. A discussion of the simulation inputs, as well as a 
program listing are included in the thesis in order to facilitate future use of the 


simulation. 
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I, INTRODUCTION 

This thesis describes an analysis of a Magnetic Anomaly Detection (MAD) 
buoy field. The analysis is based on estimates of the probability a submarine would 
be detected during an encounter with a MAD buoy. The estimates were generated 
with a program that simulates an encounter between a magnetometer based 
detection system and a submarine. The program allows a choice of two detection 
system models: a crosscorrelation detection system model and a square law detection 
system model. In the simulation, an encounter is a straight line encounter with a 
constant vertical separation between the magnetometer and the submarine. 
Probability estimates generated by the crosscorrelation detection system simulations 
are considered to represent upper bounds while probability estimates generated with 
the square law detection system simulations are considered to represent lower 
bound. 

The probability estimates generated with the simulation were used to estimate 
MAD buoy sweep widths. Probabilities of detection for a MAD buoy field were then 
estimated by using the sweep width estimates in a sonobuoy field random search 
model. Relative to operationally realizable values, probabilities estimated with the 
random search model are considered to provide lower bounds on the search 
effectiveness of a buoy field. 

Chapter 2 contains a discussion of the two detection system models, and 
Chapter 3 contains a discussion of the input parameters used in the simulation. 
Chapter 4 describes a sensitivity analysis of the simulation in which the effects of 


changes in the noise level, buoy depth and, submarine displacement are estimated. 


The computer simulation was implemented using Microsoft Quick Basic version 
4.5 on an IBM compatible 80286 computer. A program listing is contained in 
Appendix A. 














II. MODEL DESCRIPTION 


The encounter model and the computer program that was used as a basis for 


this thesis are described in Reference 1. 


A. THE MAGNETIC SIGNAL 

The encounter model is based on the following assumptions: "a submarine 
magnetic anomaly field is a magnetic dipole field that is superimposed on an earth 
magnetic field that is constant over an encounter region" [Ref. l:p. 20]. The 
magnetic signal at a sample point is determined by the Anderson formulation, and 
its duration is determined by the length of time it takes a magnetometer to pass 
through a distance that is a function of the slant range between the magnetometer 


and the magnetic dipole [Ref. 1:p. 24]. 


B. THE DETECTION MODELS 
The encounter model has two alternative binary detection models: a 
crosscorrelation detection model and a square law detection model. Both models 
are based on the following assumptions: the magnetic noise is gaussian, the noise 
samples are independent, and target course and speed remain constant during the 
detection encounter. 
1. Detection Statistics 
Let v,, v., ... » V,, represent sample values in a sample interval. Here, a 
sample value v, is the voltage generated by a magnetometer at time t,. In the 
detection models, a detector uses these sample values to determine a detection 
statistic x for the sample interval. The decision rule used by both detection models 


is: If x > x’ then the sample interval contains target signal plus noise, otherwise the 











sample interval contains noise only. Here, x’ is a value that is determined by p, the 
false alarm probability. 

For both detection models, if there is no signal present, v, = n, where n, is 
the noise present at the magnetometer and i = 1, 2,.., m. If there is signal 
present, v, = n, + 8, where s, represents the signal present at the detector. The 
average signal power S = (1/m)e3s? and the average noise power N = o” where 0 
is the standard deviation associated with the noise process. 

2. Crosscorrelation Detection Model 

The crosscorrelation detection model is based on the additional assumption 
that the detection system has knowledge of the statistical characteristics of the 
noise process and a complete prior knowledge of the input signal. This detection 
model can provide an upper bound on detection performance. It can be argued that 
a system that utilizes a computer library of signal characteristics will approach this 
upper bound. 

In this model, 

x=2v,¢ 8, 
where i = 1, 2, ... , m and the sum is over the values corresponding to a sample 
interval. When there is no signal energy present, x represents a normal random 
variable with py, = 0 and variance ox? = o2 e 3,2. The false alarm probability is 
given by: 
Pr = 1 - P(x’/ox) 

where P(x) is the standard normal cumulative distribution function. This equation 
is used to determine the value for x’ in this model. 

When the sample interval contains signal energy, x is the value of a normal 


random variable with #, = = s,2 where i = 1, 2,..., m and the sum is over the 











values corresponding to a sample interval. The result is that the probability of 
detection p, is given by: 
Pa = 1 - P(x’/oy - 2 8,2/0,). 
3. Square Law Detection Model 

The square law detection model differs from the crosscorrelation detection 
model in that the signal is determined by a stochastic process and only the 
statistical characteristics of the signal process are known. This detection model can 
provide a lower bound on detection performance. 

The statistic x for the square law detector is: 

x= v,? 
where i = 1, 2, ... ,m and the sum is over the values corresponding to the sample 
interval. 

When there is no signal present, x/o? is the value of a chi-square random 
variable with m degrees of freedom. The false alarm probability is given by: 

p, = 1- P(x’/o?|m) 
where P(x’/o?|m) is the cumulative distribution function for a chi-square random 
variable with m degrees of freedom. This equation is used to determine the value 
for x’ in this model. 

When the sample interval contains signal energy, x/o? is the value of a 
noncentral chi-square random variable with m degrees of freedom and a noncentral 
parameter = s,27/o2. The result is that for this model the probability of detection 
is given by: 

P=1- P(x'/o2|m, s,2 /o?). 
A more detailed discussion of the crosscorrelation and the square law 


detection model is provided in Reference 2. 








C. PROGRAM MODIFICATIONS 
Based on the lateral range curves that are described in Chapter II, to determine 


the probability of detecting a target transiting a MAD buoy field, the sweep width 
for each buoy was calculated. In an isotropic noise field and direct path 
environment, it is a relatively straight forward procedure to determine a sweep 
width for a sonobuoy, since the probability of detection for a submarine depends 
only on its range from the sonobuoy. However, for a MAD buoy on the other hand, 
the probability of detecting a submarine target for a MAD buoy is dependent on the 


targets course as shown in Fig. 1. 
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Figure 1. Comparisons of the Lateral Range Curves for a Target at a Constant 
Depth and Speed with Different Courses 








The program to determine lateral range curves that is described in Ref. 1 


was modified so that it would determine lateral range curves for encounters 
involving various target depths and target courses in a single computer run. In 
addition, the output of the program was altered to provide the sweep width for each 
target depth and course. 
1. Target Depth 
In the altered program, a vser can input a minimum target depth, 
maximum target depth, and an interval between calculated encounter depths. The 
depth inputs allow the user to investigate a broad spectrum of possible encounter 
depths. 
2. Target Course 
In addition, in the altered program a target’s course varies sequentially from 
000 - 150 degrees true by 30 degree increments. It is not necessary for the program 
to run the reciprocals for the courses indicated above since, for instance, the 
detection range for 000 degrees true would be the same as 180 degrees 
true. Figure 2 provides a comparison of the lateral range curve of a target on 
reciprocal headings. Notice that the two curves are mirror images of each other. For 
a target course of 180 degrees true, the lateral range curve crosses the .6 probability 
of detection line at -800 meters whereas the curve for a target course of 000 
degrees true crosses the .6 probability of detection threshold at +800 meters. 
3. Random Search Model 
A random search model described in Reference 3, pp. 18 has been 
incorporated into the computer simulation to calculate the probability of detection 
for a MAD buoy field. A random search model was chosen because other search 


models tend to be overly optimistic in their predictions [Ref 4:p. 2-1]. 
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Figure 2. Comparison of the Lateral Range Curves for a Target on Reciprocal 
Headings. 

The search model is based on the assumption that the targets motion is 

unknown and considered to be random within a defined region of area A. For a buoy 


field, the probability of detection is: 


P(det) = 1 - exp(-nfWur/A) 


where n is the number of buoys, f is the fraction of buoys that can be monitored 


simultaneously, W is the sweep width of a buoy, u is the submarine’s speed, 7 is the 


duration of the search, and A is the area of the search region [Ref. 3:p. 18]. 








A necessary condition for this formula to be justified is that the detection 


regions for the MAD buoys do not significantly overlap. If significant overlap occurs 
the actual probability of detection could be less than the calculated value. 

















III. INPUT PARAMETERS 
The input parameters for the computer program can either be entered at the 
keyboard, stored in a single file or divided into three separate files: the magnetic 
file, the signal processing file, and the kinematic file. 


A. MAGNETIC PARAMETERS 
1. The Earth Magnetic Field and Dip Angle 

Like its prototype, the modified program has the capability of either 
receiving values of the earth magnetic field H, and values of the dip angle $ as 
direct inputs or of calculating these values using an internal model. Direct input 
values for the earth magnetic field can be obtained from Defense Mapping Agency 
Hydrographic Center Chart #39 and values for the dip angle can be obtained from 
Defense Mapping Agency Center Chart #30. 

The values for the dip angle computed by the internal model will, in 
general, differ from those obtained directly from Chart #30. Table 1, which is taken 
from Reference 1, shows that a small difference in the computed and chart values 
for the dip angle does not imply that there exists a small difference between the 
computed and charted values for the earth magnetic field at the same location. For 
example, at 60°N and 30°W the difference between the charted dip angle and the 
calculated dip angle is only 1° while the difference between the charted and 
calculated earth magnetic field intensity is .12 oersted. 

2. The Submarine Magnetic Dipole 
The computer program allows the user to input permanent magnetic 


moments, target displacement, and target permeability factors or coefficients. The 
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program uses these inputs to calculate a submarines dipole moment p. Since the 


likelihood that the permanent magnetic moments for a submarine contact will be 
known is small, the values for the permanent magnetic moment are assumed to be 
zero. This approximates the condition where the submarine has been effectively 
depermed. 

The value for the target displacement used in this thesis is an average of 
displacements for various Soviet submarine classes. The information used to 


calculate the average was obtained from Jane’s Fighting Ships 1989-90 [Ref. 5]. 


Table I. SOME CORRESPONDING VALUES FOR #@ AND H.. 


Ne ee ee 
Lat. Long. Program Chart #30 Program Chart #39 


Dip ¢ Dip ¢ Field H, _—‘ Field H, 

60°N 180°W 74° 70° 63 53 
60°N 30°W 75° 74° 64 52 
30°N 90°W 62° 60° 55 51 
30°N 60°W 59° 59° 53 48 
30°N 30°E _— 7° 43° Al 42 

0° 0° -48° -25° 35 ai 
30°S 180°W -45° 55° 44 50 
30°S - 90°W_—s«-30° -34° 39 31 
30°S  90°E—s«*6;B" -66° 55 52 
60°S 30°W_—Ssé«69° 56° 60 134 


* in oersted, 1 oersted = 10° gamma 


a 


The use of this reference is not intended as a endorsement of the accuracy of the 
reference. It was used only in order to provide a convenient point of departure for 
the analysis. 

The average value d for the target displacement was calculated by summing 


the displacement d for a submarine of class i multiplied by the number n of 
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submarines in class i for all i and dividing this sum by the total number of 


operational submarines in the Soviet inventory. Using the relation: 


4 2(n;ed;) 
=n, 





the average value of a submarine contact was 5600 tons. 

The final piece of information needed is the permeability factors or the 
permeability coefficients. The permeability coefficients k, are related to the 
submarine dinnlaveriient d and the permeability factor f, by the equation: 

k, =f,ed, 
where i is the index for transverse, longitudinal, and vertical axis [Ref. 1:p. 34]. The 
computer program will accept either permeability coefficients or permeability factors 
as inputs. 

In this thesis, the permeability factors were used as an input for the 
program. The input values were obtained from Texas Instrument Report C2-61009- 
2 (Ref. 6:p. 4] and are based on a WW II vintage submarine. Care must be exercised 
when using these inputs because the units for the Texas Instrument inputs are 
foot’/ton. The program requires inputs with units of oersted-centimeter /gamma- 
ton. This problem is easily solved by dividing the Texas Instrument inputs by 3.53 
(Ref. 1:p.35]. 


B. PROCESSING PARAMETERS 

The processing parameters are crucial in that these parameters determine the 
type of magnetometer being modeled. The computer program is flexible in that it 
will accept direct inputs for the sample interval and the sample rate or it will 


calculate these inputs if they are not known beforehand. 
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1. Sampling Rate 


The program uses the value 2 e MAXF as the sampling rate where MAXF 
is defined to be the maximum magnetic signal frequency. This is the Nyquist rate 
for an ideal low pass filter [Ref. l:p. 9]. In the program, MAXF is 
2eMAXVM/MINRO where MAXVM is the maximum encounter relative speed in 
meters per second and MINRO is the minimum slant range in meters of a just 
detectable target. A just detectable target is defined to be the smallest target in 
terms of its magnetic dipole moment that can be detected with an acceptable 
probability. Since the magnetometer on a buoy is stationary, MAXVM equals the 
maximum speed of the target in knots converted to meters per second. Additional 
discussion of the sampling rate can be found in [Ref. 7:pp. 29-30). 

2. Integration Time 

The program calculates the sample interval length (integration time) as 
2eMAXRO/MINVM where MINVM is the estimated minimum encounter relative 
speed in knots converted to meters per second and MAXRO is the estimated 
maximum slant range at CPA for an encounter for the maximum expected dipole 
moment to be encountered. Additional discussion of the integration time can be 
found in [Ref. 7:pp. 28-29]. 

3. Processing Inputs 

The hypothetical MAD detection system modeled in this thesis has a 
maximum estimated slant range at CPA of 155 meters for a just detectable target. 
The maximum slant range at CPA for a target with the largest expected dipole 
moment is estimated to be 610 meters. The minimum expected encounter speed is 


5 kts and the maximum expected encounter speed is 10 kts. 
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C. KINEMATIC PARAMETERS 

The kinematic inputs include the values for magnetometer course, depth, and 
speed as well as the target course, depth, and speed. The kinematic inputs are 
basically self explanatory and with the exception of target depth, require no 
explanation. 

There are three depth parameter inputs: minimum target depth, maximum target 
depth, and depth interval. The minimum and maximum target depths define the 
range of values for which the simulation will calculate the probabilities of detection. 
The default inputs are 0 meters and 700 meters. The maximum target depth of 700 
meters was arbitrarily chosen by using Ref. 5 to determine the maximum service 
depth of current Soviet submarines. The use of 700 meters as the maximum service 
depth does not constitute a statement confirming the accuracy of the source. 

The default value for the depth interval is 50 meters. This value was 
determined to be an appropriate compromise between the longer run time that 
would result from a smaller interval and the need to provide an adequate depth 


resolution. 
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III. ANALYSIS 
The affects of false alarm rate, noise, buoy depth, and submarine displacement 


on detection are discussed in this chapter. 


A. FALSE ALARM RATE 

The simulation’s crosscorrelation and square law detection models use the input 
false alarm rate to determine the probability of false alarm p, This value represents 
the probability that a "signal present" decision will be made when in fact there is 
no signal present. As = s,2 approaches zero p, approaches p,. This relationship can 
be seen in Figure 3, at ranges greater than 1500 meters, the lateral range curves 
are converging to p,. All of the parameters for the two curves shown in Figure 3 are 
identical with the exception of the false alarm rate. In this scenario the p, generated 
by a false alarm rate of 3 false alarms per hour is .395 and the p, generated by a 
false alarm rate of .1 false alarms per hour is .013. Figure 3 illustrates the point 
that the lateral range curve with a large p, dominates the lateral range of a smaller 
P; 

For an airborne MAD detector, a false alarm rate of .3 false alarms per hour 
can be considered conservative however, for a MAD buoy field, a false alarm rate of 
3 false alarms per hour could be considered excessive while .1 false alarms per hour 
could be considered conservative. The trade off is that while a larger p, results in 
a larger p,, the higher false alarm rate leads to an increased number of false 
contacts that have to be resolved. To put this point in perspective, consider a MAD 
buoy field that consists of 10 buoys with each buoy having a life of 10 hours. A false 


alarm rate of 3 false alarms per hour for each buoy would result in, on the average, 
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Figure 3. Comparison of he Tateral Range Curves Tor Two Different False Alarm 
Rates. 






30 false alarms per hour for the buoy field. If each buoy is assumed to have a 10 


hour life then this would result in, on the average, 300 false alarms in 10 hours. On 







the other hand, a false alarm rate of .1 false alarms per hour for each buoy would, 


on the average, result in 1 false alarm per hour or 10 false alarms in 10 hours. 














B. SWEEP WIDTH 


The sweep width W of a detection system is the area under the lateral range 


curve and it is formally defined by: 


W = [pena 


where x is the lateral range and p(x) is the probability of detecting the target at a 
closest point of approach (CPA) range x during a complete straight line encounter. 
For a definite range law detection system with definite range R, the sweep width 
is W = 2eR. 

It is readily apparent that the sweep width calculated using the formal 
definition given above for the curves in Figure 3 would result in an infinitely large 
sweep width in either case. This suggests that a definite range law formula with R 
equal to the average of the .5 probability of detection ranges would be appropriate 
to define the sweep width. However, using the definite range law formula would 
conceal the asymmetry found in the magnetometer lateral range curve. For false 
alarm rates resulting in a p, of .1 or less, a value for the sweep width equal to the 
an area under the lateral range curve was used where the limits of integration 
were the ranges at which p, equals a value significantly close to p, This is a better 
choice than that above, since it can provide a more realistic representation of the 
detection process. For the analysis conducted in this thesis, the false alarm rate was 
arbitrarily set at .1 false alarms per hour which results in a p, of .013 &.:d the limits 
of integration were from the negative range to the positive range where p, equaled 


.015. Figure 4 illustrates, for a square law model, the difference in the sweep width 
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determined by the two methods of calculation. Since the lateral range curve for a 
buoy varies as a function of target course, a representative sweep width for the 
random search model was required. Based on its assumption of a uniform target 
distribution, a value was found by averaging the sweep widths that were determined 


for the target courses 000, 020, 060, 090, 120, and 150 degrees true. 
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Figure 4. Comparison Between the Sweep Width as Calculated as the Area Under 
the Lateral Range Curve vs. the 50% Probability of Detection Range. 
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C. RANDOM SEARCH MODEL SENSITIVITY ANALYSIS 


This section describes the results of a sensitivity analysis of the random search 


model. The parameters noise level, buoy depth, and target displacement were 


individually varied in the analysis for a buoy field consisting of 10 MAD buoys 
randomly placed within a 15 nm by 15 nm square. The search duration is 2 hrs and 
the target speed is 7.5kts. The fraction of the MAD buoys that can be monitored 
at one time is assumed to be 1. 
1. Noise Level 
Four different values for the magnetic noise were used in this analysis: 10 


gamma, 1 gamma, .1 gamma, and .01 gamma. As expected, an increase in magnetic 





noise lead to a significant decrease in the detection range for the buoy field, as 
shown in Figure 5. What is interesting is that the shape of the curves as a function 
of depth change as the noise levels change. 

Notice that if the magnetic noise is lowered to .01 gamma, there is not 
much variation between the buoy field probability of detections as a function of 
depth to a depth of 800 meters. This occurs because to a depth of 800 meters, the 
signal-to-noise ratio remains large and, consequently, so does the sweep width. As 
the magnetic noise increases, the variation between buoy field probability of 
detections increases due the decreasing signal-to-noise ratio. If the magnetic noise 
is increased to 1 gamma, at depths of greater than 500 meters the curve begins to 
flatten out due to the convergence of p, with p, This effect becomes even more 
noticeable when the noise is increased to 10 gamma. At depths greater than 250 
meters the buoy field detection probabilities become constant. The magnetic target 
signal has decreased to the point were, at its maximum value, p, is only slightly 


larger than p, resulting in a lateral range curve that is relatively flat. 
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2. Buoy Depth 
The buoy field probabilities of detection were calculated for three different 
buoy depths: 50 meters, 100 meters, 150 meters. Table II contains the field 
probabilities of detection for the three buoy depths. Notice that the maximum field 
detection probability occurs at the depth of the buoy and that this maximum value 
remains constant over the three depths. This result is consistent with the fact that, 
because of the magnetic signal model, only the vertical distance between the target 


and the buoy is required to determine the sweep width. 
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Figure 5. Comparison of Buoy Field Detection Probabilities with Varying Noise 
Values Using the Square Law Detection Model. 
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The final computer simulation program shown in Appendix A utilizes this 
principle to provide buoy field detection probabilities for a range of buoy depths, in 
one simulation run, without significantly increasing the simulation runtime. For 
example in Table II, if the field detection probabilities are known for the maximum 
buoy depth of 150 meters (Column 3) the field detection probabilities for all buoy 
depths less than the maximum but equal to the depth sample interval of 50 meters 
(Columns 1 and 2) are also known. 

3. Target Displacement 
This section describes a sensitivity analysis of the effects of a change in 
displacement. The maximum and minimum displacements studied were 22,400 tons 


and 700 tons. 


Table II. TABULATION CF MAD BUOY FIELD DETECTION PROBABILITIES 
FOR VARYING TARGET AND BUOY DEPTiiS. 
SE SE EEE SO PESTA 


Target Depth Buoy Depth 
meters 50 meters 100 meters 150 meters 

0 .710 .704 .697 
50 .712 .710 .704 
100 .710 .712 .710 
150 .704 .710 .712 
200 .697 .704 .710 
250 .695 .697 .704 
300 .690 695 .697 
350 .685 .690 .695 
400 .678 .685 .690 
450 .668 .678 .685 
500 .650 .668 .678 
550 .628 .650 .668 
600 .608 .628 .650 
650 591 .608 .628 
700 574 591 .608 
750 - 574 591 
800 - - 574 


Ne eg 
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Figure 6 shows buoy field detection probabilities for a square law detector. 


Notice that, over the range of operational target depths, the buoy field probability 


of detection curves for smaller displacements experience a sharper drop off in 


detection probability than the curves for larger target displacements. If target depths 


significantly greater than 700 meters had been displayed, then the same sharp 


decrease in the buoy field detection probability curves would have been observed for 


the larger displacements. 


Table III provides various statistical details about the sweep width for the 


different target displacements averaged over depths. Target depth was assumed to 
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Figure 6. Comparison of Buoy Field Detection Probabilities with Varying Target 
Displacements and Magnetic Noise of .1 Gamma. 
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Table III. STATISTICAL ANALYSIS OF THE AVERAGE SWEEP WIDTH FOR 
VARIOUS DISPLACEMENTS 


Target Displacement in tons 


22.400 11,200 5600 2,800 1,400 700 


Sample Size 15 15 15 15 15 15 
Average 2820° 2137. +1578 1106 720 446 
Median 2817 2169 1645 1185 788 516 
Mode 2804 2298 1765 1113 726 742 
Standard 

Deviation 117 154 188 212 280 280 
Minimum 2560 1796 1218 699 161 47 
Maximum 2969 2305 1774 13833 1006 763 
Interquartile 

Range 216 230 327 370 444 606 


‘All sweep widths are in meters 


be uniformly distributed between 0 meters and 700 meters. As expected the average 
of the sweep width over all depths decreases as the target displacement decreases. 
However, the range of the values over all depths increases from a minimum of 370 
meters at a target displacement of 22,400 tons to a maximum of 845m at a target 
displacement of 1400 tons. For a target displacement of 700 tons the range of values 
begins to decrease again due to the small magnetic dipole moment which results in 
the sweep width converging to a minimum value. The standard deviation follows the 
same pattern, it is inversely proportional to the target displacement. The larger the 
displacement the smaller the standard deviation. 

This relationship can be seen in Figure 6. The buoy field detection curve 
for a 22,400 ton target is relatively flat, while the curve for a 1400 ton target 
exhibits a more pronounced curve. The detection curve for a 700 ton target actually 


has an inflection point at about 500 meters target depth. This inflection point is 
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caused by the fact that the average sweep width is converging to a value greater 
than zero. 

Analysis in this section suggests that, for small submarine displacements, 
it becomes more important to place a buoy’s magnetometer near the target depth 
to ensure an even marginal probability of detection under high noise conditions. For 
example, in Figure 6, the rate of change for a 22,400 ton target (the upper curve) 
is relatively small over the operational target depths. In this case placing the buoy 
magnetometer at a depth other than the target depth does not result in a 
significant loss of detection probability. Conversely the curve for a 1400 ton target 
displays a significant rate of change and the magnetometer’s depth has a large effect 
on the ability of the buoy field to detect the target. Section four will explore the 
effects of placing the buoy magnetometers at a variety of depths within the same 
field. 

4. Multiple Buoy Depths 

In this section, the effect that placing the buoy magnetometers within the 
same field at various depths has on the detection probabilities curves is discussed. 
For this analysis the target displacement was assumed to be 1400 tons and the 10 
buoy magnetometers were placed at three different depths: three at 100 meters, 
four at 350 meters, and three at 600 meters. The size of the search area remained 
the same as in previous analysis. 

The buoy field detection probabilities were calculated using the following 
formula: 

p(det) = 2 cep, 

where ¢, is the fraction of the buoys at depth i and p, is the probability of detecting 
the target at depth i with all the buoys in the field. The fraction of the buoys at 
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the three assigned depths were as follows: .3 at 100 meters, .4 at 350 meters, and 
.3 at 600 meters. Figure 7 illustrates the results of the simulation. 

It was anticipated that placing the buoy magnetometers at a variety of 
depths within the same field would increase the probability that the buoy field 
would detect the target. It is apparent from Figure 7 that, contrary to the 
expectations, the buoy field detection probability curve for a field depth of 350 
meters dominates the multiple depth buoy field detection probability curve for all 
depths under consideration. This suggests that, assuming a uniform target depth 


distribution, deploying the buoy magnetometers to the average target depth would 
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Figure 7. Comparison of Buoy Field Detection Probability Curves for Sensor 
a of 100m, 350m, 600m, and Various Sensor Depths With In the Same Buoy 
ield. 

















result in the best probability of detecting the target. In order to confirm the results, 
a number of different scenarios were run varying buoy depth separation and target 


displacement and Figure 7 is representative of the results. 
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IV. CONCLUSIONS 


This thesis presented an analysis of a Magnetic Anomaly Detection buoy field 
that is based on data generated with a computer simulation. The simulation uses 
a crosscorrelation detection system model to provide an upper bound and a square 
law detection system model to provide a lower bound on the detection system 
performance. A random search model was used to provide a lower bound on the 
search effectiveness of the buoy field. 

Analysis has shown that the choice of the false alarm rate is the deciding 
factor in choosing the method used to calculate the sweep width. False alarm rates 
that result in a p, greater than .1 suggests that a definite range law model be used. 
For a p, that is less than or equal to .1 this assumption is not necessary since p, is 
sufficiently close to zero so that limits of integration for estimating the area under 
the lateral range can be chosen without inducing a large error. This method should 
provide a more realistic representation of the detection process. 

The random search model sensitivity analysis for noise level, buoy depth, and 
target displacement revealed that the variation between the buoy field detection 
probability as a function of depth was greater for the higher noise levels. The lowest 
noise levels showed relatively little variation in the buoy field detection probability 
over the assumed submarine operational depths, while the highest noise levels had 
a significant increase in the variation of the buoy field detection probability. The 
variability of the buoy field detection probability as a function of depth was 
decreased as the submarine displacement increased. The buoy field detection 
probability was maximized when the buoy magnetometer was at the same depth as 


the target. This indicates that the depth placement of the buoy magnetometer 
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becomes more important as noise increases or as submarine displacement decreases. 
As expected, varying the buoy depth also confirmed that the vertical separation 
between the buoy magnetometer and the target is the key factor and not the buoy 
depth. Target displacement analysis exhibited the same general trend as the noise 
level analysis. The rate of change for the large displacement targets show a less 
significant rate of change than the smaller displacement targets. Finally, analysis 
showed that varying the buoy magnetometer depth within the same buoy field does 
not improve the buoy field search effectiveness over a buoy field with all the buoy 


magnetometers at the average target depth for a uniform target depth distribution. 
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APPENDIX A. A PROGRAM LISTING 


10 CLS : CLEAR 

20 PRINT "Magnetic Anomaly Detection (MAD) Buoy Field Simulation Program" 
30 ’REM This simulation is based on a program to compute the Magnetic 

40 ’REM Anomaly Detection (MAD) lateral range function, panacea by 

50 ’REM Prof. R. N. Forrest, Naval Postgraduate School, Monterey. 

60 DIM X0(-1 TO 70), R0O(70), H(70), HS(70, 150), K(70), Puce, PDSL(71), 
AREA1T(70), AREA2T(70) 

70 PI = 3.141592654#: CON = 2 * PI / 360: KON = 1852 / 3600: N$ 
"MAD.BAS" 

80 QO = .2316419: Ql = .31938153#: Q2 = -.356563782#: Q3 = 1.781477937# 
90 Q4 = -1.821255978#: Q5 = 1.330274429# 

100 I1 = 2.515517: I2 = .802853: 13 = .010328: I4 = 1.432788: I5 = .189269: I6 
.001308 

110 PRINT 

120 INPUT "generate data or (print a program data file (g/p)"; A$ 

130 IF A$ = "G" OR A$ = "e THEN GOTO 150 

140 IF A$ = "P" OR A$ = "p" THEN GOTO 3540 ELSE GOTO 120 

150 PRINT : A$ = "a" 

160 INPUT “magnetic, processing & kinematic data entry by combined file (y/n)"; 


A$ 
170 IF A$ = "Y" OR A$ = "y" THEN GOTO 1950 
180 IF A$ = "N" OR A$ = “n" THEN GOTO 190 ELSE GOTO 160 


190 INPUT "magnetic data entry by file/keyboard (f/k)"; A$ 

200 IF A$ = "F" OR A$ = "f* THEN GOTO 810 

210 IF A$ = "K" OR A$ = "k" THEN GOTO 220 ELSE GOTO 190 
220 INPUT “latitude in decimal degrees (N +)"; LAT: LATR = LAT * CON 
230 INPUT “longitude in decimal degrees (E +)"; LNG: LNGR = -LNG * CON 
240 A "variation in decimal degrees (E +)"; DEC: DECR = DEC * CON 

250 = "a" 

260 INPUT "input dip angle (y/n)"; A$ 

270 IF A$ = "Y" OR A$ = "y" THEN GOTO 290 

280 IF A$ = "N" OR A$ = "n" THEN GOTO 320 ELSE GOTO 260 

290 INPUT "dip angle in decimal degrees (north magnetic hemisphere +)"; DIP 

300 DIPR = DIP * CON 

310 GOTO 420 

320 LATPR = 76 * CON: LNGPR = 100 * CON 

330 D = SIN(LNGR - LNGPR) * COS(LATR): E = COS(LNGR - LNGPR) * 

COS(LATR): F = SIN(LATR) 

340 X = E:Y =F 

350 GOSUB 4430 

360 T = T - (90 * CON - LATPR): E = R * SIN(T): F = R * COS(T) 

370 X = E: Y= D 

380 GOSUB 4430 


29 








390 X = F:.Y=R 
400 GOSUB 4430 
410 DIPR = ATN(2 * TAN(T)): DIP = DIPR / CON 


420 A$ = 
430 INPUT ay encounter magnetic field ey (y/n)"; A$ 
440 IF A$ A$ = "y" THEN GOTO 4 


450 IF A$ = ne OR A$ = a THEN GOTO pee ELSE GOTO 430 

460 INPUT “encounter magnetic field intensity in gamma"; HE 

470 GOTO 490 

480 HE = 70000! / pee : wa * COS(IPR) + 1) 

490 PLM = 0: PTM = 0: P 

500 A$ = “a” 

510 INPUT "input target permanent dipole moments (y/n)"; A$ 

520 IF A$ = "Y" OR A$ = "y" THEN GOTO 540 

530 IF A$ = "N" OR A$ = "n" THEN GOTO 570 ELSE GOTO 510 

540 INPUT "permanent longitudinal moment in oersted-cm3 (stern-to-bow +)"; PLM 
550 INPUT "permanent transverse moment in oersted-cm3 (port-to-starboard +)"; 
PTM 


560 INPUT "permanent vertical moment in oersted-cm3 (downward +)"; PVM 
570 INPUT "target displacement in tons"; WT 

580 A$ = "a" 

590 INPUT “input target permeability coefficients or factors (c/f)"; A$ 

600 IF A$ = "C" OR A$ = "c" THEN GOTO 620 

610 IF A$ = "F" OR A$ = "f THEN GOTO 660 ELSE GOTO 590 

620 INPUT "longitudinal permeability coefficient in cgs units"; KL 

630 INPUT "transverse permeability coefficient in cgs units"; KT 

640 INPUT "vertical permeability coefficient in cgs units"; KV 

650 FL = KL / WT: FT = KT / WT: FV = KV / WT: GOTO 710 

660 INPUT "longitudinal displacement factor in cgs units"; FL 

670 INPUT "transverse displacement factor in cgs units"; FT 

680 INPUT "vertical displacement factor in cgs units’; FV 

690 KL = ven WT: KT = FT * WT: KV = FV * WT 

700 A 

710 INPUT. "generate a magnetic data file (y/n)"; A$ 

720 IF A$ = "Y" OR A$ = "y" THEN GOTO 740 

730 IF A$ = "N" OR A$ = "n" THEN GOTO 880 ELSE GOTO 710 

740 INPUT "magnetic data file name"; M$ 

750 ON ERROR GOTO 760: GOTO 770 

760 RESUME 710 

770 OPEN “"O", #1, M$ 

tas vane #1, LAT, LNG, DEC, DIP, HE, PLM, PTM, PVM, WT, KL, KT, KV, FL, 
790 CLOSE 

800 GOTO 880 

810 INPUT "magnetic data file name"; M$ 

820 ON ERROR GOTO 830: GOTO 840 

830 RESUME 190 

840 OPEN "I", #1, M$ 

ay pee #1, LAT, LNG, DEC, DIP, HE, PLM, PTM, PVM, WT, KL, KT, KV, FL, 
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860 CLOSE 

870 DECR = DEC * CON: DIPR = DIP * CON 

880 A$ = "a" 

890 INPUT " oe data entry by file/keyboard (f/k)"; A$ 
F" OR 


900 IF A$ = " A$ = "f* THEN GOTO 1480 
910 IF A$ = "K" OR A$ = "k" THEN GOTO 920 ELSE GOTO 890 
920 A$ = "a" 


930 INPUT “input sampling period (y/n)"; A$ 

940 IF A$ = "Y" OR A$ = "y" THEN GOTO 960 

950 IF A$ = "N" OR A$ = "n" THEN GOTO 990 ELSE GOTO 930 

960 INPUT "sampling period in seconds"; DT 

970 IF DT <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT : 

GOTO 960 

980 GOTO 1120 

990 A$ = "a" 

1000 INPUT “input maximum magnetic signal frequency (y/n)"; A$ 

1010 IF A$ = "Y" OR A$ = "y" THEN GOTO 1030 

1020 IF A$ = "N" OR A$ = "n" THEN GOTO 1060 ELSE GOTO 1000 

1030 INPUT "maximum magnetic signal frequency in Hertz"; MAXF 

1040 IF MAXF <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT 

: GOTO 1030 

1050 GOTO 1110 

1060 INPUT "mir: in target slant range at CPA in meters"; MINRO 

1070 IF MINRO <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT 
: GOTO 1060 

1080 INPUT "maximum magnetometer relative speed in knots"; MAXVMK 

1090 MAXVM = MAXVMK * KON: MAXF = 2 * MAXVM / MINRO 

1100 IF MAXVM <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT 
: GOTO 1080 

1110 DT = 1 / (2 * MAXF): REM low pass filter, Nyquist sampling rate 


1130 INPUT “input integration time (y/n)"; A$ 

1140 IF A$ = "Y" OR A$ = "y" THEN GOTO 1160 

1150 IF A$ = "N" OR A$ = aa THEN GOTO 1240 ELSE GOTO 1130 

1160 INPUT "integration time in seconds"; IT 

1170 IF IT >= DT THEN GOTO 1200 

1180 PRINT : PRINT "IT = " + STR$(UT) + " seconds minimum = " + 

STR$(DT) + " seconds": PRINT 

1190 GOTO 1160 

1200 NS = 2 * INT(T / DT / 2) + 1: REM adj number of samples per integration 

time 

1210 IF NS <= 151 THEN GOTO 1370 

1220 PRINT : PRINT "IT = " + STRS$(IT) + " seconds maximum = " + 

STR$(150 * DT) + " seconds": PRINT 

1230 GOTO 1160 

1240 INPUT "maximum target slant range at CPA in meters"; MAXRO 

1250 IF MAXRO <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT 
: GOTO 1240 

1260 INPUT "minimum magnetometer relative speed in knots"; MINVMK 

1270 MINVM = MINVMK * KON 
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1280 IF MINVM <= 0 THEN PRINT : PRINT “must be greater than zero": PRINT 
: GOTO 1260 

1290 IT = 2 * MAXRO / MINVM 

1300 IF IT >= DT THEN GOTO 1330 

1310 PRINT : PRINT "IT = " + STR§$(IT) + " seconds minimum = " + 
STR$(DT) + " seconds": PRINT 

1320 GOTO 1130 

1330 NS = 2 * INTUT / DT / 2) + 1: REM adjusted number of samples per 
integration time 

1340 IF NS <= 151 THEN GOTO 1370 

1350 PRINT : PRINT "IT = " + STR$(T) + " seconds minimum = " + 
STR$(150 * DT) + " seconds": PRINT 

1360 GOTO 1130 

1370 A$ = "a" 

1380 INPUT "generate a processing data file (y/n)"; A$ 

1390 IF A$ = "Y" OR A$ = "y" THEN GOTO 1410 

1400 IF A$ = "N" OR A$ = "n" THEN GOTO 1530 ELSE GOTO 1380 
1410 INPUT "processing data file name"; P$ 

1420 ON ERROR GOTO 1430: GOTO 1440 

1430 RESUME 1380 

1440 OPEN "0", #1, P$ 

1450 WRITE #1, DT, IT, NS 

1460 CLOSE 

1470 GOTO 1540 

1480 INPUT "processing data file name"; P$ 

1490 ON ERROR GOTO 1500: GOTO 1510 

1500 RESUME 890 

1510 OPEN "I", #1, P$ 

1520 INPUT #1, DT, IT, NS 

1530 CLOSE 

1540 A$ = "a" 

1550 INPUT "kinematic data entry by file/keyboard (f/k)"; A$ 

1560 IF A$ = "F" OR A$ = "f' THEN GOTO 1770 

1570 IF A$ = "K" OR A$ = "k" THEN GOTO 1580 ELSE GOTO 1550 
1580 CM = 0: REM Magnetometer Course(CM) = 0 for a buoy 

1590 VMK = 0: REM Magnetometer Speed(VMK) = 0 for a buoy 
1600 INPUT "magnetometer depth in meters "; AM: AM = -AM 

1610 CT = 150 

1620 INPUT "target speed in knots"; VTK 

1630 INPUT "maximum target depth in meters"; AT 

1640 INPUT "minimum target depth in meters"; ATMIN 

1650 INPUT "depth interval in meters"; ATINT 

1660 A$ = "a" 

1670 INPUT "generate a kinematic data file (y/n)"; A$ 

1680 IF A$ = "Y" OR A$ = "y" THEN GOTO 1700 

1690 IF A$ = "N" OR A$ = "n" THEN GOTO 1840 ELSE GOTO 1670 
1700 INPUT "kinematic data file name"; K$ 

1710 ON ERROR GOTO 1720: GOTO 1730 

1720 RESUME 1670 

1730 OPEN "O", #1, K$ 


32 











1740 WRITE #1, CM, VMK, AM, CT, VTK, AT, ATMIN, ATINT 
1750 CLOSE 

1760 GOTO 1830 

1770 INPUT "kinematic data file name"; K$ 

1780 ON ERROR GOTO 1790: GOTO 1800 

1790 RESUME 1770 

1800 OPEN "I", #1, K$ 

1810 INPUT #1, CM, VMK, AM, CT, VTK, AT, ATMIN, ATINT 
1820 CLOSE 

1830 A$ = "a" 

1840 INPUT "generate a combined magnetic, processing & kinematic data file (y/n)"; 


A$ 

1850 IF A$ = "Y" OR A$ = "y" THEN GOTO 1870 

1860 IF A$ = "N" OR A$ = "n" THEN GOTO 2030 ELSE GOTO 1840 

1870 INPUT “combined magnetic, processing & kinematic data file name"; E$ 
1880 ON ERROR GOTO 1890: GOTO 1900 

1890 RESUME 1840 

1900 OPEN "0O"", #1, E$ 

1910 WRITE #1, LAT, LNG, DEC, DIP, HE, PLM, PTM, PVM, WT, KL, KT, KV, 
FL, FT, FV, DT, IT 

1920 WRITE #1, NS, CM, VMK, AM, CT, VTK, AT, ATMIN, ATINT, M$, P$, K$ 
1930 CLOSE 

1940 GOTO 2110 

1950 INPUT "combined magnetic, processing & kinematic data file name’; E$ 
1960 ON ERROR GOTO 1970: GOTO 1980 

1970 RESUME 1950 

1980 OPEN "I", #1, E$ 

1990 INPUT #1, LAT, LNG, DEC, DIP, HE, PLM, PTM, PVM, WT, KL, KT, KV, 
FL, FT, FV, DT, IT 

2000 INPUT #1, NS, CM, VMK, AM, CT, VTK, AT, ATMIN, ATINT, M$, P$, K$ 
2010 CLOSE 

2020 DECR = DEC * CON: DIPR = DIP * CON 

2030 INPUT ‘required false alarm rate in false alarms per hour"; FAR 

2040 PF = FAR * IT / 3600: REM false alarm probability 

2050 Y = PF: IF PF > .5 THEN Y = 1 - Y: REM inverse normal approximation 
2060 Y = SQR(LOG(1 / Y / Y)) 

2070 Y = Y- (i ee +13* Y)/(1 + Y* (4 + Y * (15 + I6 * Y))) 
2080 IF PF < .6 THEN Y = -Y 

2090 ZP = -Y 

2100 CHI = NS * (1-2/9 / NS + ZP * SQR(2 / 9 / NS)) * 3: REM inverse 
chi-square approximation 

2110 INPUT "magnetic noise in gamma"; SIG 

2120 INPUT "maximum lateral range in meters"; LRM 

2130 INPUT "lateral range step in meters"; ST 

2140 IF ST <= LRM THEN GOTO 2180 

2150 PRINT : PRINT "maximum step is " + STR$(LRM) + " meters”: PRINT : 
GOTO 2130 

2180 NC = 2 * INT(LRM / ST) + 1: REM number of lateral range function values 
2190 IF NC <= 71 THEN GOTO 2210 
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2200 PRINT : PRINT "minimum step is " + STR$(LRM / 35) + " meters": PRINT 
: GOTO 2130 

2210 ’start of mods 

2220 INPUT "number of buoys in the field"; NB 

2230 INPUT “duration of the search in hours"; SD 

2240 INPUT "search area in square nautical miles"; A: A = A * 3429904’conversion 
to square meters 

2250 INPUT "fraction of buoys that can be monitored at one time"; FR 

2260 IF FR <= 1 THEN GOTO 2280 

2270 PRINT : PRINT "fraction must be less than or equal to 1": GOTO 2250 
2280 ATZ = AT 

2290 FOR ZM = ATMIN TO ATZ STEP ATINT: REM loop for depth 

2300 AT = ZM: AREA1 = 0: AREA2 = 0 

2310 FOR ZQ = 0 TO 150 STEP 30: REM loop for course 

2320 CT = ZQ 

2330 PDCCP = 0!: PDSLP = 0! 

2340 APDCCT = 0!: APDSLT = 0! 

2350 CMR = CM * CON: CTR = CT * CON 

2360 WXK = VMK * SIN(CMR) - VTK * SIN(CTR) 

2370 WYK = VMK * COS(CMR) - VTK * COS(CTR) 


2380 Z = AM + AT: REM vertical separation (- for magnetometer below target) 
2390 X = WXK: Y = WYK: GOSUB 4430 
2400 CR = T: C = CR / CON: REM magnetometer relative course 


2410 WK = R: REM magnetometer relative speed 

2420 CMR = CM * CON: CTR = CT * CON: CR = C * CON 

2430 MCMR = (CMR - DECR): REM magnetometer magnetic course in radians 
2440 MCTR = (CTR - DECR): REM target magnetic course in radians 

2450 ILM = KL * HE * COS(DIPR) * COS(MCTR): ITM = -KT * HE * COS(DIPR) 
* SIN(MCTR) 

2460 IVM = KV * HE * SIN(DIPR) 

2470 DMX = (PLM + ILM) * SIN(MCTR) + (PTM + ITM) * COS(MCTR) 
2480 DMY = (PLM + ILM) * COS(MCTR) - (PTM + ITM) * SIN(MCTR) 
2490 DMV = PVM + IVM 

2500 X = DMX: Y = DMY 

2510 GOSUB 4430 

2520 OMLR = T: REM dipole azimuth relative to magnetic north 

2530 X = DMV: Y= R 

2540 GOSUB 4430 

2550 DM = R: OMR = T: REM dipole depression angle 

2560 OML = OMLR / CON: OM = OMR / CON 

2570 AIT = DT * NS: REM adjusted integration time 

2585 W = WK * KON: REM magnetometer relative speed in meters/second 
2590 DS = W * DT: REM distance between samples on the relative track 
2600 XO = -(NC - 1) /2* ST 

2610 FORI = 0 TONC-1 

2620 X01) = X0 

2630 X = X0: Y = Z 

2640 GOSUB 4430 

2650 RO = R: RO(I) = R: REM target slant range at CPA in meters 

2660 DELR = T: REM target depression angle complement at CPA in radians 











2670 IF RO = 0 THEN GOTO 3000: REM zero lateral range and vertical separation 
2680 DMF = DM / 10 / RO * 3: REM dipole moment factor 

2690 CMR = CR - DECR: REM target relative magnetic course 

2700 L = COS(OMR) * COS(CMR - OMLR) 

2710 M = COS(DELR) * COS(OMR) * SIN(CMR - OMLR) - SIN(DELR) * 
SIN(OMR) 

2720 N = -SIN(DELR) * COS(OMR) * SIN(CMR - OMLR) - COS(DELR) * 


SIN(OMR) 
2730 L1 = COS(DIPR) * COS(CMR) 
2740 M1 = COS(DELR) * COS(DIPR) * SIN(CMR) - SIN(DELR) * SIN(DIPR) 


2760 A2 = 2*L* L1-M * M1-N * Nl: REM Anderson Function Coefficient 
2770 Al = 3 * (N * L1 + L * N1): REM Anderson Function Coefficient 

2780 AO = 2*N* N1-L* L1-M * M1: REM Anderson Function Coefficient 
2790 SUM = 0: HMAX = 0: HMIN = 0 

2800 FOR J = 0 TONS - 1 

2810 S = (J - (NS - 1) / 2) * DS: BA = S / RO: REM Anderson Function Argument 
2820 AF = 1 / (1 + BA * BA) * 2.5: REM Anderson Function Factor 

2830 HSF = (A2 * BA * BA + Al * BA + AO) * AF: REM magnetic signal factor 
2840 HS(I, J) = DMF * HSF: REM magnetic signal value 

2850 IF HS, J) > HMAX THEN HMAX = HSI, J) 

2860 IF HS(, J) < HMIN THEN HMIN = HS, J) 

2870 SUM = SUM + HSF * HSF 

2880 NEXT J 

2890 H(I) = HMAX - HMIN 

2900 K(I) = SQR(SUM) 

2910 VV = -ZP + DMF * SQR(SUM) / SIG 

2920 LAM = DMF * DMF * SUM / (SIG * SIG): AA = NS + LAM: BB = 1 + LAM 
/ (NS + LAM) 

2930 ZN = -SQR(2 * CHI / BB) + SQR(2 * AA / BB - 1): X1 = VV 

2940 GOSUB 4490 

2950 IF Y1 > 1 THEN Y1 = 1 

2960 PDCC(I) = Y1: X1 = ZN 

2970 GOSUB 4490 

2980 IF Y1 > 1 THEN Y1 = 1 

2990 PDSL(I) = Y1 

3000 XO = XO + ST 

3010 REM this section calculates the sweep width as the area under 

3020 ’the lateral range curve 

3030 IF I > 0 THEN 

3040 IF PDCCP <= PDCC(I) THEN 

3050  AREAI1 = (PDCCP * ST) + (ST * ABS(PDCC(1) - PDCCP)) / 2) 


(PDCC(D) * ST) + (ST * ABS(PDCCP - PDCC(I))) / 2) 
3080 END IF 


3090 REM APDCCT & APDSLT are the areas under their respective lateral 
3100 ’range curves 

3110 APDCCT = APDCCT + AREAI1 

3120 IF PDSLP <= PDSL(I) THEN 

3130 AREAI2 = (PDSLP * ST) + (ST * ABS(PDSL() - PDSLP)) / 2) 


2750 N1 = -SIN(DELR) * COS(DIPR) * SIN(CMR) - COS(DELR) * SIN(DIPR) 


: 
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3140 ELSE 

3150 AREAI2 = (PDSL() * ST) + (ST * ABS(PDSLP - PDSL())) / 2) 

3160 END IF 

3170 APDSLT = APDSLT + AREAI2 

3180 END IF 

3190 PDCCP = PDCC(1) 

3200 PDSLP = PDSL(D 

3210 NEXT I 

3220 C = C / 360: C = (C - INT(C)) * 360 

3230 IF C < 0 THEN C = 360 + C 

3240 OML = (OML + DEC) / 360: OML = (OML - INT(OML)) * 360 

3250 IF OML < v THEN OML = 360 + OML 

3260 REM AREA1&2 sums the sweep width over all courses for a given depth 

3270 AREA1 = AREA1 + APDCCT: AREA2 = AREA2 + APDSLT 

3280 NEXT.ZQ 

3290 AREA1T(COUNT) = AREAI / 6 calculating the average sweep width 

3300 REM PDET1TQ( calculates the buoy field probability of detection 

3310 for a crosscorrelation detector and PDET2T() for a square law detector 

3320 PDET1T(COUNT) = 1 - EXP(-NB * FR * AREALT(COUNT) * VTK * 1852 * 

SD / A) 

3330 AREA2T(COUNT) = AREA2 / 6 

3340 PDET2T(COUNT) = 1 - EXP(-NB * FR * AREA2T(COUNT) * VTK * 1852 * 

SD / A) 

3350 COUNT = COUNT + 1 

3360 NEXT ZM 

3370 PRINT : A$ = "a" 

3390 INPUT "generate a program data file (y/n)"; A$ 

3400 IF A$ = "Y" OR A$ = "y" THEN GOTO 3420 

3410 IF A$ = "N" OR A$ = "n" THEN GOTO 3650 ELSE GOTO 3390 

3420 INPUT "program data file name"; D$ 

3430 ON ERROR GOTO 3440: GOTO 3450 

3440 RESUME 3420 

3450 OPEN "0", #1, D$ 

3460 WRITE #1, LAT, LNG, DEC, DIP, HE, PLM, PTM, PVM, WT, KL, KT, KV, 

FL, FT, FV, DT, IT, AIT 

ay are #1, NS, CM, VMK, AM, CT, VTK, AT, ATMIN, ATINT, COUNT, FAR, 
, ,S 

3480 WRITE #1, NB, SD, A, FR, LRM, DS, NC, M$, P$, K$, E$ 

3490 FOR I = 0 TO COUNT 

3500 WRITE #1, AREAIT(I), AREA2T(I) 

3510 NEXT I 

3520 CLOSE 

3530 GOTO 3650 

3540 INPUT "program data file name"; D$ 

3550 ON ERROR GOTO 3560: GOTO 3570 

3560 RESUME 3540 

3570 OPEN "I", #1, D$ 

3580 INPUT #1, LAT, LNG, DEC, DIP, HE, PLM, PTM, PVM, WT, KL, KT, KV, 

FL, FT, FV, DT, IT, AIT 
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3590 INPUT #1, NS, CM, VMK, AM, CT, VTK, AT, ATMIM, ATINT, COUNT, FAR, 
PF, SIG, ST 

3600 INPUT #1, NB, SD, A, FR, LRM, DS, NC, M$, P$, K$, E$ 

3610 FOR I = 0 TO COUNT 

3620 INPUT #1, AREAIT(I), AREA2T(I) 

3630 NEXT I 

3640 CLOSE 

3650 PRINT : A$ = "a" 

3660 INPUT "print encounter parameter values (y/n)"; A$ 

3670 IF A$ = "Y" OR A$ = "y" THEN GOTO 3690 

3680 IF A$ = "N" OR A$ = "n” THEN GOTO 4140 ELSE GOTO 3660 
3690 LPRINT 


3700 LPRINT “program file name "+ N$ 
3710 LPRINT “program data file name "+ D$ 
3720 LPRINT "magnetic data file name "+ M$ 
3730 LPRINT "processing data file name "+ P$ 

3740 LPRINT "kinematic data file name "+ K$ 
3750 LPRINT "combined magnetic, processing & kinematic data file name " + E$ 
3760 LPRINT "encounter latitude (decimal degrees) "» SPC(2); LAT 
3770 LPRINT "encounter longitude (decimal degrees) "; SPC(2); LNG 
3780 LPRINT “encounter variation (decimal degrees) "; SPC(2); DEC 
3790 LPRINT "encounter dip angle (decimal degrees) "; SPC(2); DIP 
3800 LPRINT "encounter magnetic field intensity (gamma) ", SPC(2); HE 
3810 LPRINT "permanent longitudinal moment (oersted-cm3) ", SPC(2), 
PLM 

3820 LPRINT "permanent transverse moment (oersted-cm3) ", SPC(2); 
PTM 

sae LPRINT "permanent vertical moment (oersted-cm3) ", SPC(2); 
3840 LPRINT "target displacement (tons) ", SPC(2); WT 
3850 LPRINT "target longitudinal permeability coefficient "; SPC(2); KL 
3860 LPRINT "target transverse permeability coefficient ", SPC(2); KT 
3870 LPRINT "target vertical permeability coefficient " SPC(2); KV 
3880 LPRINT "target longitudinal permeability factor "» SPC(2); FL 
3890 LPRINT "target transverse permeability factor ", SPC(2); FT 
3900 LPRINT "target vertical permeability factor "; SPC(2); FV 
3910 LPRINT "sampling period (seconds) *» SPC(2); DT 
3920 LPRINT "integration time (seconds) " SPC(2); IT 
3930 LPRINT "adjusted integration time (seconds) ", SPC(2); AIT 
3940 LPRINT "number of samples per encounter "» SPC(2); NS 
3950 LPRINT "magnetometer depth (meters) ", SPC(2); AM 
3960 LPRINT "target speed (knots) "» SPC(2); VTK 
3970 LPRINT "maximum target depth (meters) ", SPC(2); AT 
3980 LPRINT "minimum target depth (meters) "; SPC(2); 
ATMIN 

3990 LPRINT "depth interval (meters) " SPC(2); ATINT 
4000 LPRINT "distance between samples on the relative track (meters)"; SPC(2), DS 
4010 LPRINT "false alarm rate (false alarms per hour) "; SPC(2); FAR 
4020 LPRINT "false alarm probability ", SPC(2); PF 
4030 LPRINT "magnetic noise (gamma) "» SPC(2); SIG 
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4040 LPRINT "maximum lateral range (meters) ", SPC(2); LRM 


4050 LPRINT “lateral range step (meters) " SPC(2); ST 
4060 LPRINT “number of lateral range function values " SPC(2); NC 
4070 LPRINT “number of buoys in the field "" SPC(2); NB 
4080 LPRINT “duration of the search (hours) "» SPC(2); SD 
4090 LPRINT "search area in square nuatical miles ", SPC(2); A 


4100 LPRINT "fraction of buoys that can be monitored at one time °; SPC(2); FR 
4110 FOR I = 0 TO 21 

4120 LPRINT 

4130 NEXT I 

4140 PRINT : A$ = 

4150 INPUT ery ie field probability of detection values (y/n)"; A$ 
4160 IF A$ = "Y" OR A$ = ‘y" THEN GOTO 4180 

4170 IF A$ = "N" OR A$ = "n" THEN GOTO 4260 ELSE GOTO 4150 
4180 LPRINT D$; " buoy field probability of detection values" 

4190 LPRINT : LPRINT 

4200 LPRINT "depth p(cc) p(sl)" 

4210 LPRINT "meters" 

4220 LPRINT 

4230 FOR I = 0 TO COUNT - 1 

4240 LPRINT ATINT * I; TAB(10); PDET1T(); TAB(24); PDET2T(1) 
4250 NEXT I 

4260 PRINT : A$ = "a" 

4270 INPUT "print average sweep width values (y/n)"; A$ 

4280 IF A$ = "Y" OR A$ = "y" THEN GOTO 4300 

4290 IF A$ = "N" OR A$ = "n" THEN GOTO 4380 ELSE GOTO 4270 
4300 LPRINT D$; " average sweep width values" 

4310 LPRINT : LPRINT 


4320 LPRINT “depth avg(cc) avg(s})" 
4330 LPRINT "meters meters meters” 
4340 LPRINT 


4350 FOR I = 0 TO COUNT - 

4360 LPRINT ATINT * I; TAB(10), AREAIT(I); TAB(24); AREA2T(I) 
4370 NEXT I 

4380 PRINT : A$ = 


4390 INPUT mcontinde-t the program (y/n)"; A$ 

4400 IF Ag = "V" OR Ag = * " THEN GOTO 10 

4410 IF A$ = "N" OR A$ = "n" THEN GOTO 4420 ELSE GOTO 4390 

4420 END 

4430 R = SQR(X * X + Y * Y): REM rectangular to polar conversion 

4440 IF R = 0 THEN T = 0 0: RETURN 

4450 IF apa = 1 THEN @ = SGN() * (PI / 2) ELSE Q = ATN(XX / R / 
SQR(1-X*X/R /R 

1460 IF ABS / By) « ) THEN T = (PI / 2) * (1 - SGN(Y)) ELSE T = (PI / 2) - 
ATN(Y / R / SQR(1- Y * Y / KR / R)) 


4470 IF Q < 0 THEN T = 2* PI-T 

4480 RETURN 

4490 Y1 = X1: IF X1 < 0 THEN Y1 = -Y1: REM normal approximation 
4500 G = 1/ (1 + QO * Yl) 
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4510 Yl = EXP(-Y1 * Y1 / 2) / SQR(2 * PI) *G * (Q1 + G*(Q2+G*(Q3+G 
* (Q4 + G * Q5)))) 

4520 IF X1 >= 0 THEN Yl = 1- Yl 

4530 RETURN 
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